# -*- coding: utf-8 -*-
"""
Created on Tue Jul 20 22:41:52 2021

@author: LSY
"""

from __future__ import division
import numpy as np
import scipy.io as spio

# ======================================================================================= #
# ==================================  Sub 1 ============================================= #
                                 # CSTR 1 and 2
# ======================================================================================= #
                                 
  
                               
def ode_sub1_4cstr(x_scale,u,w): # contains Reactor 1 Reactor 5 and Material 1
    
    xs_sub1 = [2.78, 363, 2.58, 356, 2.6, 355, 2.6, 392]
    Nx_sub = np.array([8])
    min_x_sub1 = np.zeros(Nx_sub[0]) 
    max_x_sub1 = np.zeros(Nx_sub[0])  
    lb = 0.5
    ub = 1.5
    for i in range(Nx_sub[0]):
        min_x_sub1[i] = lb*xs_sub1[i]
        max_x_sub1[i] = ub*xs_sub1[i]   
    delta_x_sub1 = max_x_sub1 - min_x_sub1  

    
    x = x_scale*delta_x_sub1 + min_x_sub1  
    
    C01 = 4;     C02 = 2;       C03 = 3;       C04 = 3.5;
    F_01 = 5;    F_02 = 10;     F_03 = 8;      F_04 = 12;
    V_1 = 1;     V_2 = 3;       V_3 = 4;       V_4 = 6;
    F1=35;       F2=45;         F3=33;         Fr1 = 20;   Fr2 = 10;   
    E_1=5e4;     E_2=7.5e4;     E_3=7.53e4;    R=8.314;
    T_01=300;    T_02=300;      T_03=300;      T_04=300;
    D_H1=-5e4;   D_H2=-5.2e4;   D_H3=-5e4;
    k_10=3e6;    k_20=3e5;      k_30=3e5;
    C_p=0.231;   rho=1000;

    alpha_1=D_H1*k_10/rho/C_p;
    alpha_2=D_H2*k_20/rho/C_p;
    alpha_3=D_H3*k_30/rho/C_p;
    
#   States and parameters for sub 1
    x1 = x[0];      x2 = x[1];     x3 = x[2];     x4 = x[3];     # CA1, T1, CA2, T2 
    x5 = x[4];      x6 = x[5];     x7 = x[6];     x8 = x[7];    
            
#   Inputs
    Q1 = u[8];     Q2= u[9];     Q3 = u[10];    Q4 = u[11];

#   Dynamics
    dx1=F_01/V_1*(C01-x1)+Fr1/V_1*(x3-x1)+Fr2/V_1*(x7-x1)-(k_10*np.exp(-E_1/R/x2)+k_20*np.exp(-E_2/R/x2)+k_30*np.exp(-E_3/R/x2))*x1
    dx2=F_01/V_1*(T_01-x2)+Fr1/V_1*(x4-x2)+Fr2/V_1*(x8-x2)-(alpha_1*np.exp(-E_1/R/x2)+alpha_2*np.exp(-E_2/R/x2)+alpha_3*np.exp(-E_3/R/x2))*x1+Q1/rho/C_p/V_1
    dx3=F1/V_2*(x1-x3)+F_02/V_2*(C02-x3)-(k_10*np.exp(-E_1/R/x4)+k_20*np.exp(-E_2/R/x4)+k_30*np.exp(-E_3/R/x4))*x3
    dx4=F1/V_2*(x2-x4)+F_02/V_2*(T_02-x4)-(alpha_1*np.exp(-E_1/R/x4)+alpha_2*np.exp(-E_2/R/x4)+alpha_3*np.exp(-E_3/R/x4))*x3+Q2/rho/C_p/V_2
    dx5=(F2-Fr1)/V_3*(x3-x5)+F_03/V_3*(C03-x5)-(k_10*np.exp(-E_1/R/x6)+k_20*np.exp(-E_2/R/x6)+k_30*np.exp(-E_3/R/x6))*x5;
    dx6=(F2-Fr1)/V_3*(x4-x6)+F_03/V_3*(T_03-x6)-(alpha_1*np.exp(-E_1/R/x6)+alpha_2*np.exp(-E_2/R/x6)+alpha_3*np.exp(-E_3/R/x6))*x5+Q3/rho/C_p/V_3;
    dx7=F3/V_4*(x5-x7)+F_04/V_4*(C04-x7)-(k_10*np.exp(-E_1/R/x8)+k_20*np.exp(-E_2/R/x8)+k_30*np.exp(-E_3/R/x8))*x7;
    dx8=F3/V_4*(x6-x8)+F_04/V_4*(T_04-x8)-(alpha_1*np.exp(-E_1/R/x8)+alpha_2*np.exp(-E_2/R/x8)+alpha_3*np.exp(-E_3/R/x8))*x7+Q4/rho/C_p/V_4;    

    
#   extend states  # scale (1./1800)*1./delta_x_sub1*
    dxdt_sub1 =1./delta_x_sub1* [
                               dx1 + w[0],
                               dx2 + w[1],
                               dx3 + w[2],
                               dx4 + w[3],
                               dx5 + w[4],
                               dx6 + w[5],
                               dx7 + w[6],
                               dx8 + w[7],
                              ]
    return np.array(dxdt_sub1) 
#
            
# ======================================================================================= #

def measurement_sub_1(x):  #这里要不要进行归一化？
    
    xs_sub1 = [2.78, 363, 2.58, 356, 2.6, 355, 2.6, 392]    
#    xs_sub1 = [2.78882185, 363.41439929,  2.58905237,   356.54409766,   5, 10, 1, 3, 10]  
    Nx_sub = np.array([8])
    min_x_sub1 = np.zeros(Nx_sub[0]) 
    max_x_sub1 = np.zeros(Nx_sub[0])  
    lb = 0.5
    ub = 1.5
    for i in range(Nx_sub[0]):
        min_x_sub1[i] = lb*xs_sub1[i]
        max_x_sub1[i] = ub*xs_sub1[i]
    
    delta_x_sub1 = max_x_sub1 - min_x_sub1 
    
    C_sub1 = np.array([[1,0,0,0,0,0,0,0],
                       [0,0,0,0,0,0,0,1]])     
    
#    C_sub1 = np.array([[0,0,0,0,0,0,0,1],
#                       [0,0,0,0,0,0,1,0],
#                       [1,0,0,0,0,0,0,0]])
    
    diag_delta_x_sub1 = np.diag(delta_x_sub1) # make the scales in a diagonal matrix
    C_sub1_scaled = np.dot(C_sub1,diag_delta_x_sub1)
#    C_matrix = np.identity(145)
    y_sub1 = np.dot(C_sub1_scaled,x) + np.dot(C_sub1,min_x_sub1)

    
#    y_sub1 = x[0]
    
    return y_sub1  